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A non-ideal MHD Gadget: Simulating massive galaxy clusters 
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ABSTRACT 

Magnetic fields in the intra-cluster medium of galaxy clusters have been studied in the past 
years through different methods. So far, our understanding of the origin of these magnetic 
fields, as well as their role in the process of structure formation and their interplay with the 
other constituents of the intra-cluster medium is still limited. In the next years the up-coming 
generation of radio telescopes is going to provide new data that have the potential of setting 
constraints on the properties of magnetic fields in galaxy clusters. 

Here we present zoomed-in simulations for a set of massive galaxy clusters (M^ ^ 
\Q^^h~^MQ). This is an ideal sample to study the evolution of magnetic field during the 
process of structure formation in detail. Turbulent motions of the gas within the ICM will 
manifest themselves in a macroscopic magnetic resistivity ry™, which has to be taken explic- 
itly into account, especially at scales below the resolution limit. We have adapted the MHD 
GADGET code by Dolag & Stasyszyn (2009) to include the treatment of the magnetic re- 
sistivity and for the first time we have included non-ideal MHD equations to better follow 
the evolution of the magnetic field within galaxy clusters. We investigate which value of the 
magnetic resistivity rjm is required to match the magnetic field profile derived from radio ob- 
servations. We find that a value of ?],„ ^ 6 x 10^^ cm^s^^ is necessary to recover the shape 
of the magnetic field profile inferred from radio observations of the Coma cluster This value 
agrees well with the expected level of turbulent motions within the ICM at our resolution limit. 
The magnetic field profiles of the simulated clusters can be fitted by a /3— model like profile 
(Cavaliere & Fusco-Femiano 1976), with small dispersion of the parameters. We find also 
that that the temperature, density and entropy profiles of the clusters depend on the magnetic 
resistivity constant, having flatter profiles in the inner regions when the magnetic resistivity 
increases. 

Key words: (magnetohydrodynamics)MHD - magnetic fields - methods: numerical - galax- 
ies; clusters 



1 INTRODUCTION 

Magnetic fields are an important ingredient to understand the phys- 
ical processes taking places in the intra-cluster medium (ICM) of 
galaxy clusters. Their presence is demonstrated by radio obser- 
vations, which, since the last 30 years, have revealed diffuse and 
faint radio sources fill ing the central Mpc'^ of some galaxy cluster s 
(radio halos, see e.g. iGiovannini et al.ll2009l : IVenturi et alJuOOSh . 
These sources arise because of the interaction of highly relativis- 
tic electrons with the ICM magnetic fields. About 30 radio halos 
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are known so far, and all of them are found in clust ers with clear 
signatures of on-going or recent merge r activity (e.g. iBuotdllOOll : 
lOovoni et alj200lMCassano et alj201(ll) . The origin of the relativis- 
tic particles still needs to be understood, although several models 
have been proposed. Shocks and turbulence associated with merger 
events are expected to inject a considerable amount of energy in 
the ICM, that could compress and amplify the magnetic field and 
(re-)accelerate rela tivistic electrons, giving thus rise to the observed 
radio emission (see lFerrari et al.l2008l . lDolag et al . 2008, for recent 
reviews of the subject). Understanding the magnetic field amplifi- 
cation and evolution during the process of structure formation is 
then mandatory for modeling the acceleration, transport and inter- 
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actions of non-thermal energetic particles and thus to understand 
the observed emission. In addition, an accurate modeling of the 
magnetic field properties is necessary to understand both the heat 
transport and the dissipative processes in the ICM. 
The properties of magnetic fields in the ICM have been inves- 
tigated in the past through cosm ological sim ulations, perf ormed 
with different numerical codes JDolag et alJ |l999. 20 021, l2005l: 
Dubois & Tevssiej|2008l . lOolag & Stasvszvnll2009l . ICollins et alj 



(Ml 



20IC) and also 
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served data is necessary to constrain the main magnetic field prop- 
erties, and it is starting now to be feasible thanks to the progress 
that has been done in the recent years. One key aspect is that, so 
far, large scale radio emission is mainly detected in very massive 
clusters. Such massive systems are not easily studied by numerical 
simulations, since the size of the density fluctuations responsible 
for the formation of massive halos is large, i.e.~ 20 /i~^Mpc, and 
the value of the cosmological parameter ag in the standard AC DM 
model requires that statistically a total volume of ~ 200h~^ Mpc^ 
needs to be sampled by simulations in order to produce at least one 
cluster as massive as ~ lO^^ft^^A/0. An important step for study- 
ing non-thermal phenomena is to perform simulations based on ex- 
tremely large cosmological volumes, e.g. I Gpc side-length. Such 
large volumes cannot be simulated at the resolution reached by ob- 
servations, so that re-simulation tech niques ha ve been developed 
(e.^.G RAFIC Bertschingei 199j; ZIC lTormen et al. 1997; Je nkinsl 
I2OICI) . When such high resolution is reached, the physic of the 
baryonic component must be followed with particular care. The 
magnetic field amplification, in particular, depends on the small 
scale motions of the gas. Hence, as the resolution increases smaller 
scale motions are revealed, a nd the magnetic fi eld amplification in- 
creases accordingly (see e.ie. lDolag et alj2008r) . 
In this paper we present a set of galaxy clusters extracted by a low 
resolution DM simulation and re-simulated at high resolution (the 
softening length is ~5 kpc h^^) within a cosmological framework 
in order to resolve scales comparable to those reached by observa- 
tions. This work is focused on the 24 most massive galaxy clus- 
ters (A/200 > 10^^h~^Mpc) of our sample. Simulations are per- 
formed for the first time relaxing the assumption of ideal MHD, 
and including a resistivity term in the induction equation (rim). Our 
sample of simulated galaxy cluster is publically available for fur- 
ther studieo In this paper we present the simulated cluster sample: 
the MHD implementation with some test problems (Section. [2l(, the 
re-simulation technique used (Section. [3] and more detailed in the 
appendix); the effect of different values for rim are analyzed and 
discussed in Section [4] where the main properties of the clusters 
are also presented. Finally, discussion and conclusions are reported 
in Section[5] 

This is a first paper aimed at presenting the cluster sample, the 
zoom-in initial conditions, and the non-ideal MHD implementa- 
tion in the GADGE T code. This sample has also been used by 
iFabian et al.l j201lh for a study of the scaling relations of X-ray 
mass proxies. In a future paper the authors will investigate in more 
detail the cluster properties, and the interplay between thermal and 
non-thermal components in the ICM. 



2 NON-IDEAL MHD SIMULATIONS 

Within the last decade, cosmological simulations of structure for- 
mation have shown that the observed properties of magnetic fields 
in galaxy clusters are direct consequences of turbulent amplifi- 
cation driven by the the s tructure formation process ( Dolag et alj 
J999, 2002, 2005; Brugge n et ai]|2005l; JPubois & Tevssieill2008l ; 
Dolag & Stasvszvn, ,200^ ; ICollins et alj l2010h . Simulations per- 



formed with different codes reach good agreement in predict- 
ing that the ratio of the bulk kinetic energy to the thermal en- 
ergy has an upper limit of ~ 10-20% (see e.g. the review by 
iBorgani & Kravtsovll2009l and references therein). R ecently, non- 
cosmological MHD simulations of merging galaxies ( Kotarba et al.l 
2009, 2010) predict that the magnetic field is amplified up to a level 
close to ~ 10-20% of the thermal energy. The same is expected for 
the ICM of galaxy clusters. Although the properties of magnetic 
fields in galaxy clusters are still not strongly constrained from the 
observational point of view, present data suggest that the magnetic 
field energy content is not amplified up to the level of the kinetic 
energy. In the Coma cluster, for e xample, the turb ulent energy con- 
tent is ~10% of the thermal one l lSchuecker et al...2004. ). whereas 
the magnetic energy c ontent associated with in the observed mag- 
netic field of 4.7/iG teonafed e et"ai]|2010 l) is only 1.6% of the 
thermal one. Dissipative processes could possibly explain the satu- 
ration of magnetic fields far below the level of the kinetic energy. 
Such dissipative processes, driven by the physical properties of the 
ICM plasma, are no t investigate in numerical simulations so far, but 
iDolag & Stasvszvnl ([2009) have shown that dissipative processes 
driven by numerical diffusivity may alter the central properties of 
the magnetic field profiles obtained by numerical simulations. 
The simulations we present in this paper were carried out 
with GADGET-3 JSpringell |2005|) , the current version of the 
parallel TreePM-l-SPH simulation code GADGET JSpringel et al.1 
[2001'). It includes an entropy-conserving formulation of SPH 
(ISpringel & Hernquist 2002), the implementation of ideal MHD 
( iDolag & S tasvszvn 2 009|) and a n implem entation of a diverg ence 
cleaning scheme (Dol ag & Stasv szvn 2OO9;! B0rve et alj200lh . The 
cosmological simulations presented here assume an initially homo- 
geneous magnetic field of lO^^'^G co-moving. 
In previous works it was usually assumed that the electric conduc- 
tivity of the gas is infinite, meaning that the second term of the 
induction equation (Eq.[Tll vanishes (rjm ~ 0). 
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= VX(UXB) + VX (77mV X B). 



(1) 



This assumption results in a magnetic field frozen into the gas. We 
have extended the treatment of the induction equation to cover the 
resistive MHD equation. Here we will assume for simplicity a spa- 
tially constant resistivity term rim. In Section l4~2l the physical ori- 
gin of this term is analyzed and the assumption will be discussed. 
Under the constraint V • -B = 0, and rim spatially constant, the 
induction equation for resistive MHD can then be written as: 



dB 

'dt 



{B-V)v~B{\/ ■v)+ri„ 



(2) 
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The resistivity dependent terms have been implemented in the 
co de follow ing the approach ad opted for the artificial dissipation 
by Price & Monaghan (20Q4a,b. 120050 . In particular, we refer to 
iDoIag & Stasvszvr> ([2009) where the artificial dissipation term has 
been implemented in the GADGET code. More specifically, the re- 
sistive term is included in the induction equation as 
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Figure 1. Comparison of the results from the simulations (diamonds) to the 
analytic solution (lines) at different output times. The magnetic resistivity 
r\jyi was set to 1 in this test. For graphical reasons, only one diamond each 
8th particle in x-direction was plotted for every time step. 
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Where i and j refer to two generic particles in the simulation, ri,j 
is the distance between particle i and j, W is the SPH kernel, and 
the factor (Ha^)^^ = J| takes into account that the internal time 
variable in GADGET is the expansion parameter a. The resistivity 
term implemented in the induction equation causes a change in the 
Entropy A at the rate 
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where Wij indicates the mean between the two kernels Wi 
and Wj, 7 is the adiabatic index of the gas. We refer to 
iDolag & Stasvszvnl f 2009) for more details about the numerical im- 
plementation. Since we basically replaced the artificial dissipation 
already implemented and tested iDolag & Stasvszvn 2009 ) with a 
physically motivated magnetic resistivity, the only tests that are 
left to be performed to validate the numerical scheme are those 
regarding the ability of the code to reproduce the correct dissipa- 
tion timescale. This can be done, in the case of a spatially constant 
Tim, by investigating the magnetic field evolution for simple test 
problems. 



2.1 Test 1: A one-dimensional slab in a 3D setup 

We consider first the time evolution of a one-dimensional mag- 
netic field (B = B(t)y) in a one dimensional slab at rest having 
side length 4L. In order to test the code within the configuration 
used for cosmological simulations, we performed the test in a 3D 
setup using a glass-like particle di stribution and solving a planar 
test problem within this 3D setup ( IWhitdll996h . We started with 
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Figure 2. Comparison of the results from the simulations (diamonds) to the 
analytic solution (lines) at different output times. The magnetic resistivity 
rjm was set to 1 in this test. For graphical reasons 1 diamonds each 4th 
particle in x-direction is plotted for every time. 



700 X 10 X 10 particles, having a mean inter-particle separation 
along the x axis of 5.7xlO~^L. Using 64 neighbors within the 
SPH interpolant this correspond to roughly 35 resolution elements 
per length L. 

The induction equation here reduces to 

A^B 



= r)r, 



dB 
'dt 



which has the analytical solution: 



(5) 



B{t) = exp -T]mt 



/27r\2\ „ /23;7r 



(6) 



We set Bx = B and By = Bz = 0, and followed the evolution of 
an initial magnetic field B = _Bosin(27ra::/(4L)). 
In Figure [T] the time evolution of the system is shown. The results 
obtained from the numerical simulation (diamonds) are compared 
with the analytic solution (lines) for various time steps, showing an 
excellent agreement. 

2.2 Test 2: Magnetic diffusion across a step in a 3D setup 

As a second test, we consider here a one dimensional slab. The 
magnetic field is described by _B = B{x, t)y, and a step profile for 
the magnetic field was included according to: 



B{x,0) 



+Bo x>0 
-Bo 2; < 



(7) 



As in the previous test, the simulation was performed in a full, three 
dimensional setup using a glass-like particle distribution and solv- 
ing a planar test problem within this 3D setup. We started with 
700x10x10 particles, having a mean inter-particle separation along 
the X axis of 5.7 x 10""^!/. Using 64 neighbors within the SPH inter- 
polant this correspond to roughly 35 resolution elements per length 
L. Under the constrain 



B{~L,t) = -B{L,t) = Bo 



(8) 
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meaning that the magnetic field is held fixed at two points (±1/) the 
solution of the diffusion equation can be written as 



Table 1. Properties of the high mass cluster set. 



B(x,t)^Bo — A > —exp 
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(see IWilmot-Smith et alJuOOSl) . In Figure |2] the results of the nu- 
merical simulation (diamonds) are compared to the analytic solu- 
tion (lines) at different time steps, as reported in the figure panel. 
The magnetic field diffuses rapidly and converges towards the 
steady-state solution, B{x) = Bo{x)/L. Since we have not imple- 
mented the necessary boundary conditions to keep B — Bo fixed 
at the borders, this simulation was stopped early. 

2.3 Total energy conservation 

The two tests described in Section. 12.11 and 12.21 demonstrate 
the ability of the code to correctly solve the diffusion equation. 
When real physical problems are considered, the magnetic energy 
dissipated is explicitely added to the energy equation, i.e. it is 
transfeiTed to the system, so that the total energy is conserved. 
Hence, within the cosmological simulations presented in the 
following Sections, the energy of the dissipated magnetic field 
is transferred into heat. This energy is added explicitly to the 
internal energy, similarly to what is done when artificial magnetic 
resistivity is used as a regularization scheme to suppres s numer- 
ical instabilities ( Price & Monaghan 2007, Dolag & Stasvszvnl 
1200^1) . We refer to Dolag & StasvszviJ J2009i) for the details 
of the numerical implementation, and in p articular, to section 
3.1 an d Figure 4 (upper middle panel) of iDolag & Stasvszvnl 
( 1200 9) where such issues are discussed and analyzed. In the 
two tests presented in Section 12.11 and 12.21 the conversion of 
the dissipated magnetic energy into heat has been switched off, 
since the solution we compare with do not include such conversion. 



3 CONSTRUCTING THE CLUSTER SET 

3.1 The parent simulation 

The clusters were selected from a N-body cosmological simulation 
performed according to a flat ACDM cosmological model, with 
f2m =0.24 (the matter density parameter), Qtat- =0.04 (the contri- 
bution given by baryons), h =0.72, and as = 0.8. The power spec- 
trum for the primordial density fluctuations P{k) oc k" is char- 
acterized by n — 0.96. This simulation was carried out with the 
massively parallel T REE+SPH code GADGET- 3, the new version 
of the GADGET code JSpringel et al. 2001; Springel 2005) and con- 
sists of a periodic box of size 1 Gpc h^^ . The cluster identification 
was p erformed at z = using a standard Friend of friends algo- 
rithm JDavis et alJI 19851) . The linking length was fixed to 0.17 the 
mean inter-particle separation between DM particles, correspond- 
ing to the virial over-density in the adopted cosmological model. 
This large simulated cosmological box contains 64 clusters with 
Mfof > lO'^^ h~^MQa.l z = Q. Hence, it represents a proper 
sample to study the general properties of massive galaxy clusters. 
Since we want to analyze the magnetic field properties, and com- 
pare our results with those found from Coma cluster observations, 
a statistical set of galaxy clusters with masses similar to the one 
of Coma is needed. Note that up to now, no such sample of high 
resolution re-simulations of massive galaxy clusters has been con- 
structed. 
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Col. I: Cluster name; Col. 2: Total mass inside -Rviri 

Col. 3: Virial radius; 

Col 4: Estimated X-ray Luminosity in the band 0.1- 10 keV; 

Col 5: Mean temperature (mass weighted); 

All quantities are coinputed inside R^i^- 

3.2 Cluster selection and Initial Conditions 

Clusters were selected from the parent simulation on the basis of 
their mass only. We selected the 24 most massive objects among 
those with A/fof > 10^^ ft" '^M© and re-simulated each of these 
clusters at highe r resolution by usin g the Zoomed Initial Condi- 
tions code (ZIC, iTormen et al.ll 19971) . In the appendix the iterative 
procedure used to obtain the high resolution initial conditions is de- 
scribed in detail. The setup of initial conditions was optimized to 
guarantee a spherical volume around each cluster with radius of ~ 
5-6 virial radii (.Rvir) simulated at high resolution (HR region) and 
free of contamination by low resolution, boundary particles. Two 
of the cluster initially selected by the Friend of Friends algorithm 
turned out to have a companion with mass > 10^^ h^^AlQ. Other 
systems are undergoing a merger event at 2 = with less mas- 
sive companions. In addition, other clusters with masses between 
10" /i"^Af0and 10^^ /i"^MQwhere found in the HR region of 
the main targets. 50 of them are cleaned by low resolution parti- 
cles inside their virial radius. Therefore, the final sample consists of 
76 clusters with masses larger than IO^'^H'^Mq, comprising both 
isolated and merging systems. The massive cluster set is shown 
in Table 12.31 In the Appendix (Table lAlt more details about the 
cluster surroundings are given, while in Table I A2] the clusters with 
M > lO^*/i~^M0, found within 5 Rvir from the massive targets, 
are listed. They are simulated at high resolution up to 1-5 Rvh, as 
reported in that Table. 

The virial mass of each cluster was defined as the mass contained 
within a radius encompassing an average density equal to the virial 
density, pvir, predicted by the top-hat spherical collapse model. For 
the assumed cosmology it is p^ir = 95pc, where pc is the critical 
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Figure 3. Projected X-ray surface brightness of the clusters in our sample computed in the range 0.1- 10 keV (square root scale). The side of each box 
corresponds to ~ 1.6 X 1.4 iJvir- 



cosmic density JEke et al J 19960 . In this work we focus onto the 24 
originally selected clusters, as they represent a statistical well de- 
fined, volume and limited sample of massive clusters. 
Once the ICs for the DM components have been obtained, gas par- 
ticles were added (see appendix for details). The mass of DM and 
gas particle is O.84x^/i~^Af0and O.16xlO'-'/i'^A/0respectively. 
The gravitational softening length used is 5 kpc h^^, which corre- 
sponds to the smallest SPH smoothing length reached in the dense 
cluster centers. 



3.3 The high mass cluster set 

Simulations of these cluster set including radi ative losses and 
star-formation are presented in lFabian et al.l ( l201ll) . Here, we focus 
on non-radiative simulations, since our aim is to study the effect of 
the magnetic field, and of non-ideal MHD. 

From the final snapshots of these simulations we derived the 
projected X-ray surface bright ness images, by using a map-making 
algorithm (IDolag et alj l2005h . The predicted emission of every 
SPH particle is projected along the line of sight considering an 



6 A. Bonafede, K. Dolag, F. Stasyszyn, G. Murante, S. Borgani 



integration depth of ± 5 i?vir around the center of simulated 
clusters. The X-ray Luminosity {Lx) is computed in the range 
0.1-10 keV. The X-ray surface brightness images of the clusters are 
shown in Figure [3] The values of Lx and of the gas temperature 
inside the virial radius are reported in Table 12.31 Clusters in 
different dynamical state belong to this sample and consequently 
the X-ray surface brightness images show quite different mor- 
phologies. Several clusters are disturbed in the very internal 
part, indicating that a merger event has just occurred (e.g.D_12), 
while other clusters have multiple peaks in the X-ray images, like 
e.g.DA. Some clusters appear to have a regular shape, and others 
are going to interact with a smaller halo, that is visible in the 
X-ray images (e.g.DA). In the sample we also found an ongoing 
merger event between two massive clusters (D_8, interacting 
with another cluster of M > 10^^/i"^Mo). We note that the 
over all range of morphologies found in this mass limited sample 
compares qualitatively well with complete, obs erved samples (like 
the REXCESS sample, iBohringer et al.l 120070 . where also such 
extremely perturbed systems are found). 
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Table 2. Diffusion coefficients used in our simulations (in internal and phys- 
ical units) together example values of turbulent length-scales and velocities 
which would correspond to such values. 



3-5 xlO^s 
3x10^5 
1 X 1029 
2x10^'' 
2xl02« 
6x102^.5x10^9 



Process ref. 

CR propagation (value at 1 GeV) 1 

CR driven dynamo in galaxies 2 

powering the Coma radio halo 3 

turbulent cascade observed in Coma at 2 kpc 4 

turbulent cascade from simulations at 7.8 kpc 5 

iron abundance profile in clusters 6 



Table 3. Values of the diffusion coefficients common ly used in th e litera - 



ture. and observationally inferred. References are: 
2: Lesch & Hanasz (2003), 3: Schlickeiser et al. IT98' 
l2004b . 5: llVIaieretalJ<2009l) . 6: lRebusco et alj<2006l 



Strong et al. ('2007^ 
, 4:|Schuecker et aj 
See text for details. 



4 MAGNETIC FIELDS IN MASSIVE CLUSTERS 

The properties of the ICM magnetic fields start now to be bet- 
ter understood, thanks to an increasing effort in analyzing Fara- 
day Rotation Images of sour ces located either i nside clusters 
and in their background (e .g Murgiaet all 12004 , Iciarke 20041 



Johnston-HoUitt et al. 



2006LlBonafede et al.l 



nd (e. gjfVLurgia et atj UUU4). lUiarKei |2UU41 
20041 IVogt & EnBlinI |2005| iGovoni et al.l 



20IO t). In general, the magnetic field in clus- 



ters inferred from these observations is found to be consistent with 
a magnetic field driven by the turbulence within the ICM and gen- 
erally shows a radial decline. Once the density profile p{r) has 
been inferred from X-ray observations, the magnetic field profile 
in galaxy clusters is supposed to follow the gas density profile ac- 
cording to: 



B(r) = i3op(r)". 



(10) 



The fluctuations within the magnetic field are usually modeled as- 
suming a power-law power spectrum, described by a slope rj, a 
maximum length scale Amax (which can be related to the outer 
scale of the turbulence in within the ICM) and a minimum length 
scale Amin (which in case it is resolved, could be related to dissi- 
pative scales, either viscous or resistive). These model parameter 
are inferred by comparing the expected Rotation Measure statis- 
tics (mean, dispersion, auto-correlation function and structure func- 
tion) and the polarization properties of the radio galaxies to the ob- 
served ones. So far the magnetic field in the Coma cluster is the 
one that is best constrained. It has been inferred from RM obser- 
vations of seven radio-sources located at projected distances of 50 
to 1500 kpc from the cluster center. The best fit model results to 
be the one with Bo = 4,.7toltJ.G, a = 0.51°;?. and Amin ~ 2 
kpc JBonafede et al.l2OI0r) . Although previous cosmological MHD 
simulations of galaxy clusters produced magnetic field configura- 
tion which le ad to Rotation Measure s tatisti cs similar to the ob- 
served ones (Dola g etalj|l999ll2002ll2005l), the m agnetic field 
profile tended to be steeper, with a ~ 1 dDolag e t al. 2001). In 
addition, the values of the central magnetic field obtained from 
high-re solution simulations resulted to be slightly larger than ob- 
served dDonnert et al.ll2009f) . but it was noticed that the magnetic 
field profiles are significantly altered if the underlying numerical 



MHD i mplementation suffers fro m the presence of numerical dif- 
fusion dDolag & Stasvszvij|2009l) . 



4.1 Testing the effect of the magnetic resistivity 

Having a stable numerical scheme at hand, which does not suffer 
from numerical diffusion outside the SPH smoothing length (Sta- 
syszyn 2011, in preparation), we can investigate for the first time 
the role of a physically motivated resistivity -qm in shaping the ICM 
magnetic field profile. From our set of massive clusters, which have 
all masses comparable to the Coma's one, we selected four ob- 
jects that at jz = show X-ray morphologies similar to the one of 
Coma. In particular, we avoid selecting clusters with very spherical 
morphology as well as clusters with clear multiple X-ray bright- 
ness peaks. Figure|4] shows the X-ray morphology of the 4 selected 
clusters for the 3 spatial projection directions. This sub-set of clus- 
ters has been simulated with different value of rim, with the aim 
of studying the resulting shape and central value of the magnetic 
field. Figure |4] shows the magnetic field profiles of those clusters 
compared to the best fit model for the Coma cluster, encompassed 
by the ± 3cr region. Whereas all magnetic field profiles obtained 
from the simulations are within the 3(7 region in the outer parts, the 
profiles with small magnetic diffusion (rim =1.5x10^^001^8"^) 
are always above this region towards the center. For large mag- 
netic diffusion (rim = 6x lO^^cm^s"^) half of the simulated pro- 
files are above the best fit model, the other half below the best fit 
model in the central part. From that, we conclude that a value of 
Vm — 6x lO^^cm^s"^ (20 in the code internal units) is the one that 
provides the best match with that inferred from Coma cluster obser- 
vations. The numerical diffusion inside the SPH smoothing length 
is of the order of lO^^cm^s"^ that is several orders of magnitudes 
lower than the one we have implemented as magnetic resistivity, 
and thus does not affect our results. 



4.2 Physical origin of the magnetic resistivity 

In the previous Sections we have shown that a relatively large value 
of rim is required in the induction equation (Eq.|2} to match the ra- 
dial profile and the central value of the magnetic field inferred from 
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Figure 4. Top: The right panel shows the Coma cluster X-ray surface brightness from the ROSAT All Sky Survey, in the energy band 0.1- 2.4 keV (color 
coded). The shown region size is ~ 3 X 3 Mpc, corresp onding to ~ 1 i?vir x 1 ^vir- The left panel shows the magnetic field profile that gives the best fit to 
Faraday rotation observations (see lBonafede et alj|2010l for details). Bottom: Projected X-ray suiface brightness for the clusters used to test the value of r]m.- 
(Top: D_2 (left), D_13 (right); Bottom: D_5 (left), D_20 (right))x. X Y and Z projection are shown in the upper right, lower left and right sub-panels of each 
panel. The size of the projected X-ray images corresponds to 2 iJvir- The magnetic field radial profile ({B)) is plotted in the upper left panel for values of 
r]m =1.3, 3 and 6 X 10 cm s^^ (5, 10 and 20 in code units). Magenta lines indicate the profile of the Coma cluster as inferred from radio RM observations. 
The shaded line repre sents the best fit, the dotted-shaded lines indicate the flatter and shallower profile that are compatible within 3cr with radio RM data (see 
iBonafede et alj2010l for further details). 
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Coma cluster observations. In order to correctly interpret this result, 
it must be kept in mind that the induction equation (Eq.|2} describes 
the evolution of a magnetic field B at our resolution limit (which 
is of order of 10 kpc). The turbulent cascade is expected to develop 
down to smaller scales, where unresolved turbulent motions would 
contribute to the diffusion described by rim ■ Hence, we can define 
the diffusion coefficient rim as 



^m — ?7Coulomb + ^^turb, 

with rycouiomb related to the thermal conductivity a by 

'/Couloinb . 

47rcr 



(11) 



(12) 



Following ISpitzen J1956r) . when the mean free path is determined 
by Coulomb collisions, the thermal conductivity of the plasma can 
be expressed as 



2 1/2 

7re m. 



(47reo)2(feT)3/2 



In(A), 



(13) 



A being the Coulomb logarithm. For a typical cluster environment 
(e.g.densities n « 10~^cm~'' and temperatures T ~ 10**°K), the 
diffusion coefficient is: 



^culomb 



2 X io"r 



-3/2 



(14) 



Hence, the diffusion coefficient, arising from the gas thermal con- 
ductivity, does not significantly contribute to the evolution of B in 
the induction equation. On the other hand, in a turbulent plasma the 
motion of charges will be a random walk char acterized by a length 
scale A turb and by a velocity uturb- Following l Oennis & ChandranI 
( 120051) . the plasma turbulent diffusion coefficient ?7turb can be de- 
fined as: 



?7turb ~ 0.1 Aturb X Uturb 



(15) 



Typical values of vturb at our resolution of several kpc, correspond- 
ing to scales Xturh that fall below our resolution limit, will be sev- 
eral tens of km s~^, and will lead to diffusion coefficients simi- 
lar to the one that we used in our simulations (see Table |2j. Es- 
timates of Vturb at these small spatial scales cannot be provided 
by any observation so far. However, it is possible to infer such es- 
timates from the values of v obtained at larger spatial scales, as- 
suming that the turbulent power spectrum can be described by a 
singl e power-law down to th e small scales of interest. Using X-ray 
data, ISchuecker et alj ( l2004f) derived pseudo-pressure fluctuations 
maps of the gas in the Coma cluster. They revealed the presence of a 
scale-invariant pressure fluctuation spectrum, that is consistent with 
the Kolmogorov slope, and could estimate the size of the turbulent 
eddies in the range from 40 kpc to 100 kpc. On smaller scales, 
the number of photons detected were not sufficient for a reliable 
pressure measurement. The energy content associated with these 
turb ulent motions is estima ted to be roughly 10% of the thermal 
one JSchuecker et alj|2004h . The sound velocity within the Coma 
cluster (T ~ lO^K) is ~ 1500 kms"\ Therefore the turbulent ve- 
locities associat e d with the largest scales (~ 100 kpc) found by 
ISchuecker et al.l ( 120041) would correspond to ^ 470 km s^^ . As- 
suming a Kolmogorov-like power spectrum, this translates into a 
turbulent velocity of ~ 30 kms~^ at a length scale of 2 kpc, that 
is the minimum scale revealed by Rotation Measure observations 
([Bonafede et al. 2010). A turbulent velocity of ~ 30 kms~^ at 2 
kpc would yield to rim ~ 2 x lO^^cm'^s"^, similar to the value we 
have used in the simulations. A sample of clusters for which Vturb, 
Xturb, and the power spectrum slope are estimated observationally 
would of course allow us a better and more reliable comparison. 



Although such observations are not available in t he literature so 
far, an other estimate for riturb has been derived bv lRebusco et al.l 
( 120061) . The authors have analyzed the effect of turbulent diffusion 
on the iron abundance profiles in the ICM for a sample of clusters, 
finding riturb in the range 6 x 10^^ - 4.5 x 10^^. 
Estimates of Vturb and Xturb can also be deriv ed by c osmological 
simulation. Diffe r ent numerical schemes (e.g. lOolag et al. 200: 
IVazza et al.ll20o3 : llapichino & Niemeveill2008l : IVazza et siBw' 
have been optimized to follow the evolution of tur bulent flows 
withi n the ICM of simulated galaxy clusters (see also IZhuravleval 
uOllh . These works indicate that the energy in turbulent motions is 
~ 10 — 20 % of the thermal one at 2: = within the virial radius. 
In particular, iMaier et al.l ( 120091) have measured the spectral prop- 
erties of the gas velocity field, finding a good agreement with the 
Kolmogorov power spectrum slope over scales ranging from 300 
kpc down to the scale correspondent to the Nyquist frequency. The 
velocity of the turbulent eddies at scales of 10 kpc is estimated to 
be ~ 50-100 km s~^, resulting is rim ~ x 10^*, in good agreement 
with t he values adopted i n the simulations presented here. In ad- 



2 X 10^** cm^s 



dition, iMaier et al.l ( 120091) have found ryturb 

using Aturb ~ 7.8 kpc h~^ and fturb ~ 60 kms~^. It is also 
worth mentioning that the value of the diffusion coefficient rim is 
within the range 7?^ »; 3 x lO^^cm^s"^ -i^m ~ 3x 10^^ cm^s"^ 
which are the values nee ded to operate a cosm ic ray driven dy- 
namo within a galaxy (see lSieikowski et al.l2O10l) and to power the 
Coma radio halo by an in-situ acceleration model respectiveljlf] (see 
iSchlickeiser et al. 1987). 

All these different values for rim are reported in Table |3] In sum- 
mary we can conclude that our inferred value of rim ~ 6 x 
lO^^cm^s"^ at our unresolved scales of ~ 10 kpc is well in range 
with what would be expected from turbulent motions within the 
ICM. 



4.2.1 About the use of a constant rim 

The Equations \TT\ and [TS] clarify the physical origin of the resis- 
tivity term that we have implemented in the induction equation. In 
this work, the value of rim has been kept constant throughout the 
whole re-simulations. It is clear that having a rim that changes as a 
function of Vturb and Xturb locally would allow one to follow the 
evolution of the magnetic field more properly. Identifying the turbu- 
lent motions to compute the most correct value of rim at every step 
during the simulation is however not feasible. Different algorithms 
have bee n developed to identify and analyze the turbulent motion s 
(see e.g. iDolag et alj|2005l : iMaier et all [20091: IVazza et aklboill) . 
These algorithms need to subtract large-scale laminar motions be- 
fore revealing the turbulent patterns, and can then be applied in the 
post processing once the simulation is run. We can verify the va- 
lidity of the assumed constant value of rim by checking which val- 
ues of Aturb and vturb are obtained at different distances from the 
cluster center by the above mentioned works. ( IMaier et alj|2009l) 
have computed the profile of the turbulent velocity for a simulated 
galaxy cluster. The velocity profile, once scaled at the length scale 
of 7.8 kpc - the highest resolved region of the simulation- shows 
a rather flat profile, with values ranging from ~50 to ^^90 Kms~^ 
within the cluster virial radius. This implies that the assumption of 



^ We note however that the re-acceleration model proposed in that 
work is due to Alt'en modes, while more recent works indicate that 
the re-acceleration is due to MHD modes at very small scales, see 
e.g lBrunetti & Lazariai] J201lh . 



Cluster name 


Bo 


re 
kpc 


fJ- 


x' 


D.1 


4.7 


339 


0.40 


0.9 


BJ. 


6.8 


295 


0.52 


0.7 


DJ 


2.5 


361 


0.33 


1.8 


D.4 


2.7 


285 


0.43 


1.2 


D^ 


3.0 


346 


0.35 


1.7 


D.6 


5.1 


414 


0.57 


0.9 


D.7 


6.5 


362 


0.58 


1.1 


D.8 


3.5 


342 


0.34 


1.0 


DS 


5.4 


404 


0.49 


0.9 


D.IO 


3.9 


319 


0.38 


1.9 


D_ll 


3.9 


293 


0.31 


2.0 


D_12 


3.3 


415 


0.51 


1.4 


D.13 


6.5 


332 


0.59 


1.0 


D.14 


4.6 


329 


0.44 


0.7 


D_15 


2.9 


344 


0.30 


1.9 


D_16 


5.6 


354 


0.54 


0.9 


D.17 


5.5 


431 


0.43 


1.1 


D.18 


6.4 


327 


0.63 


1.1 


D_19 


2.6 


341 


0.29 


1.7 


D_20 


4.3 


323 


0.63 


0.7 


DJ21 


2.8 


311 


0.43 


1.9 


DJ22 


4.6 


423 


0.46 


0.9 


D_23 


9.9 


241 


0.53 


1.2 


DJ24 


4.9 


375 


0.55 


1.0 


Mean values 


4.7±1.7 


346±47 


0.46±0.11 




Coma Cluster 


4 7+0.7 
^■'-0.8 


291 ±17 


n ,r;+0.17 





Table 4. Results of the fit of the magnetic field profiles to Eg. 1161 



a constant r]m, although not optimal, is well justified in our case. 
The simulations presented here represent a good starting point to 
investigate for the first time the effects related to such resistivity 
term. 



4.3 Magnetic properties of massive clusters 



We finally simulated the whole cluster sample, fixing the magnetic 
resistivity to our inferred value of rjm = 6 x lO^^cm^s^^. There- 
fore, we can study, for the first time, the scatter of the magnetic 
field properties in massive clusters using a volume limited sample. 
Figure [5] shows the mean magnetic field within 0.3 x Rvii, cor- 
responding to roughly 1 Mpc for our set of massive clusters. The 
simulations scatter mildly around the value inferred from observa- 
tions of Coma. Our mean value is 2.6/iG with and rms of 0.6/iG. 
We point out here again that the sample of massive clusters com- 
prises objects that have very different dynamical state al z — 0. It 
is interesting to note that the mean magnetic field, averaged over 
the central Mpc"^, does not depend on the present dynamical state 
of the cluster at z — (see also Section|5]) 



4. 3. 1 Magnetic field profiles of massive galaxy clusters 

In Figure[6]the magnetic field profiles are shown for all the clusters 
in the sample. In the right panel of that Figure the mean and the dis- 
persion of the magnetic field profiles are compared with the best fit 
for the Coma cluster. It is worth stressing that the exact shape of the 



profile inferred from Coma observations is given ad hoc as a para- 
metric model to fit the data. Hence, it is not clear how significant 
the differences between simulations and observations in the exact 
shape are. Nonetheless, the fit to the observations lies completely 
within the scatter of the profiles predicted by our simulations. This 
is a non-trivial result, confirming previous findings that the mag- 
netic field in galaxy clusters is shaped by the (turbulent) motions 
within the ICM and therefore reflects a natural prediction of the 
structure formation process. 

Although the mean magnetic field profile shows a good agreement 
with the one inferred from Coma observations, there are differ- 
ences in the shape of the individual profiles, likely reflecting the 
dynamical state and the different morphologies of the individual 
objects. Magnetic field profiles are usually compared to the gas 
density profiles, in order to derive a scaling with the radial dis- 
tance from the cluster center. Here we adopt another approach and 
fit the magnetic field profiles dire ctly to a "/3 model-like" profile 
JCavaliere & Fusco-Femiandll976h . that for magnetic fields is usu- 
ally written as: 



(16) 



B(r) = Bo[l + 



where Bo, re, and ji are free parameters. The fits have been per- 
formed in the range 50-2000 kpc, to properly compare with the 
results obtained from Coma observations. The results of the fit are 
shown in Table|4] The mean values of Bo and /x are reported is the 
last row of the Table and compared to those of the Coma cluster. 
While previously numerical simulations indicated a steeper profile 
of the magnetic field with the gas density (i.e.a ~ 1 in Eg.llOt. 
and thus with the radial distance from the cluster center, now the 
effect of the magnetic resistivity is that of flattening the profiles, 
reaching a better agreement with observations. In the case of the 
Coma cluster, a value of a = O.SIqi gives the best fit with the 
observations, corresponding to /i ~ O.SS^q Jg, in very good agree- 
ment with the mean of the best fit for our simulated clusters, that is 
fj, ~ 0.46 ± 0.11. Hence, not only the mean value of B over the 
central Mpc^ has a small dispersion in this high-mass cluster set, 
but also the central value of Bo, and its slope with the gas density, 
as derived from the beta-model fit, are quite similar. 



4.4 Magnetic field and thermal properties 

We present in this Section a first overview of the thermal properties 
of the ICM in the presence of a diffusive magnetic field. A more 
detailed analysis will be performed in a second paper, where the 
whole sample will be analyzed. Here we investigate if and how the 
presence of resistive magnetic field may affect the ICM properties 
of the four clusters we have simulated with different values of rjm ■ 
In Figure |7] the density, temperature, and entropy profiles of these 
clusters are shown for different values of rim ■ The profiles converge 
at distances larger than few % of the virial radius, while differ in the 
very inner region of the clusters. As mentioned in Section 1231 the 
magnetic field energy that is dissipated during the cluster formation 
is transformed into heat. Although the dynamical effect of a mag- 
netic field of the order of ~ I- 10 fiG in the cluster cores is negligi- 
ble, the overall effect of the magnetic force and pressure integrated 
over a Hubble time results in a change of the density and tempera- 
ture profile. As the resistivity constant rjm increases, the amount of 
magnetic field energy, that is dissipated and hence converted into 
heat, increases accordingly. An additional source of heating is then 
present in the cluster central region, that has the effect of flatten- 
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Figure 5. Left: Magnetic field averaged over the central 0.3 Rvir versus virial mass of the our cluster set (Blu e diamonds). The red cross refers to the mean 
magnetic field for the Coma cluster The error-bar refers to the 3a" of the chi^ given bv lBonafede et al.l 1201011 . Right: Magnetic field in the cluster center, as 
results from the fit of a /3-model profile, versus the cluster virial mass. 
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Figure 6. Left: magnetic field strength profile for the whole cluster set. Right: in blue mean (continuous lin e) and dispersion (dot- dashed) of the magnetic field 
strength profile. Red dashed line refers to the best fit obtained from RM observations for the Coma cluster jBonafede et alJlOlOl) . 



ing the temperature profile. The higher pressure that would result 
from a higher temperature is then balanced by reducing the gas den- 
sity in the cluster central region, up to a factor 2. The temperature 
and density profiles do not change adiabatically, as demonstrated 
by the entropy profiles of the clusters. The entropy, computed as 



fcT/i 



,2/3 



flattens in the inner region of the clusters, indicat- 



ing that the transport of low entropy gas in inhibited. 
In Figure[8] the magnetic field within 0.3 x 7?vir, is plotted versus 
the cluster mean temperature computed over the same region. Al- 
though the sample is small, and the value of the magnetic fields 
within 0.3 x i?vir varies of a factor ~ 2, a trend is suggested. 
Magnetic field in higher temperature clusters seem to be higher. 
The correlation should be better investigated with a higher sam- 
ple of simulated galaxy clusters, since observational data do not 
sugg est a trend of the R M in clusters depending on the tempera- 
ture JGovoni et alJuOlOl) . We note also that such a trend is much 



less visible when the value of Bo, resulting from the /?-model fit is 
compared with the cluster mean temperature (FigurejS] right panel). 



5 DISCUSSION AND CONCLUSIONS 

We have presented a set of simulated galaxy clusters. It consists 
of 24 massive objects (Afvir > lO^^/i~^M0) re-simulated at 
high resolution up to 5-6 virial radii, plus 50 more clusters with 
i\fvir > lO^''/i~^Af0that fell within this high resolution region. 
This large set permits to study the cluster properties in a wide range 
of masses and at high resolution (see e.g. Fabjan et al. 2011). The 
evolution of the clusters has been follo wed using the MHD imple - 
mentation within the GAD GET- 3 code JDolag & StasYszvnll2009h . 
that has been here modified in order to include the magnetic resis- 
tivity term in the induction equations. It is the first time that this 
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Figure 7. Density (left column), temperature (middle column) and entropy (right column) profiles for the clusters D_2, D_5, D_13 and D_20 from top to 
bottom respectively. Different colors refer to different values of the magnetic resistivity constant Tjm, as indicated in the panels. 



term is analyzed in the context of cluster formation and evolution. 
In this first paper we have presented the zoomed initial conditions, 
the non-ideal MHD implementation, and the average properties of 
the more massive clusters when a magnetic resistivity term is in- 
cluded in the MHD equations. Further analysis will be performed 
in a future paper (Paper II, Bonafede et al. in prep) where the phys- 
ical implications will be discussed in more detail. 
Our main results can be summarized as follows: 



• The magnetic field profiles obtained with non-ideal MHD can 
reproduce the profile inferred from Faraday Rotation Measures of 
the Coma cluster. Four clusters having X-ray morphologies similar 
to the one of the Coma cluster have been selected to test the effect 
of changing the constant rim used in the induction equation. The 
best agreement with the limits given by observations is achieved 
with 77™ = 6 X 10^'' cm^s"\ 



• Non-ideal MHD equations have been implemented within the 
GADGET code. The tests performed on two different problems show 
that the numerical implementation is accurate and can be used to 
study the effect of the magnetic resitivity. 



• The whole sample has been simulated using rjm — 6 x 
lO^^cm^s"^, and the derived magnetic field profiles are consis- 
tent with the Coma profile. The best-fit found for the Coma profile 
lies in fact between the rms of the simulated profiles. 
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Figure 8. Left: Magnetic field averaged over the central 0.3 iJvir versus the mean temperature of the cluster set (Blu e diamonds). The red cross refers to 
the mean magnetic field for the Coma cluster. The error-bar refers to the 3cr of the chP given bv lBonafede et al.l I2010h . Right: Magnetic field in the cluster 
center, a s results from the fit of a /3-model profile versus the cluster mean temperature inside 0.3 i?vir . The temperature for the Coma cluster is the one 
given by lAmaud et alj 1200 ih . computed inside 0.25 of the Coma virial radius. All the quantities are computed from the dissipative simulation runs with 
)7^ = 6 X lO^'^cm^s-l . 



• We have fitted the magnetic field profile with a /3-like model, 
finding that the magnetic field profile of the simulated galaxy clus- 
ters can be well reproduced by values Bo — 4.7 ± 1.7, and 
fi = 0.46 ± 0.11 (see Ea.ll6t. in good agreement with the value 
found for the Coma cluster. The value of /i would correspond to a 
value of Q ~ 0.6 (Eg.llOt for a Coma-like gas density profile. 

• We have investigated possible correlations of the magnetic 
field strength with the cluster mass. The magnetic field strength, av- 
eraged over a central spherical volume of O.SRvirh'^ in radius, is 
similar for all the cluste rs in the sample, i n agre ement with what has 
been recently found bv lBonafede et al.l <201lh . This indicates that 
the presence of radio halo emission, found in a fraction of massive 
galaxy clusters, cannot be attributed to a difference in the magnetic 
field strength. A mild dependence of the magnetic field strength 
with cluster temperature is indicated by these simulations. 

• The density, temperature and entropy profiles of the simulated 
galaxy clusters have been derived for different values of rim ■ We 
find that the effect of a magnetic diffusive constant is visible in 
such profiles, leading to flatter temperature and entropy profiles in 
the inner region of the cluster {R ^ 0.1-Rvir at maximum). 

The cluster sample and the new MHD-implementation we have pre- 
sented is suitable to investigate other issues that are not discussed 
here, and that will be studied in a future paper, such as the interplay 
of the magnetic field with the thermal gas of the ICM (e.g. how is 
the thermal conduction modified, the role of the magnetic pressure 
in suppressing the cooling in the inner regions). In the next years, 
the LOw Frequency ARray (LOFAR) and the Expanded Very Large 
Array (EVLA) will allow us to improve our knowledge of the non- 
thermal component of the ICM, and a larger sample of data will be 
soon available for a more complete comparison. 
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APPENDIX A: GENERATING THE ZOOMED INITIAL 
CONDITIONS 

Al Dark matter re-simulations 

Creating zoomed initial conditions is essential to extend the dy- 
namical range accessible through cosmological simulations, which 
is needed to study the detailed structure of objects, like e.^. galaxy 
clusters, with appropriate resolution. Since hydrodynamic simula- 
tions are sensitive to boundary conditions, regions around galaxy 
clusters have to be re-simulated with high resolution as well. In the 
last years the peripheral regions around galaxy clusters are also at- 
tracting more and more interest, given the increased sensitivity of 
modem instruments. Here we have optimized our initial conditions 
to study a statistical sample of massive clusters with reasonable 
computational resources. Our procedure is based on the ZIC code 
JTormen et al.lll997h and we describe here the iterative procedure 
that we have used to obtain such highly optimized, zoomed initial 
conditions for our cluster sample. 

We started from a large, cosmological, dark-matter only simula- 
tions, performed according to the 'concordance' ACDM cosmo- 
logical model (Oa = 0.76, fio =0.24, h =0.72 and an = 0.8). The 
spectral index of power spectrum for the primordial density fluc- 
tuations (P{k) oc fc") is n = 0.96. This simulation, that we refer 
to as 'parent simulation', was carried out with the massively paral- 
lel TREE-l-SPH code GADGET-2 (Springel 2005) and consists of 
a periodic box of 1 h^^ Gpc size. The cluster identification was 
performed at 2: = using a standard Friend of Friends algorithm. 
The linking length was fixed to 0.17 of the mean inter-particle sep- 
aration between DM particles, to reflect the virial over-density for 
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Table Al. 



Cluster i? cleaned J^dm N of nearby clusters 

[flvir] withM> 1014]/i-1Mq 



D.l 


5.2 


1.618 


D.2 


5.4 


1.518 


D.3 


5.3 


1.49 


D.4 


5.4 


1.482 


D^ 


5.0 


1.537 


D.6 


5.0 


1.165 


D.7 


5.4 


1.776 


D.8 


5.3 


1.993,1.170 


D.9 


5.2 


1.657 


D.IO 


5.1 


1.705 


D.ll 


5.3 


3.163 


D.12 


5.5 


1.678 


D.13 


5.6 


1.171 


D.14 


6.0 


1.557 


D.15 


5.5 


1.840 


D.16 


5.2 


1.385 


D.17 


5.5 


1.813 


D.l 8 


5.1 


1.356 


D.19 


5.1 


1.316 


D.20 


5.2 


1.067 


D.21 


5.9 


1.623,1.011 


D.22 


5.2 


1.674 


D.23 


5.1 


1.880 


D.24 


5.0 


1.507 



Col. 1 : Cluster name; Col. 2: Number of virial radii cleaned by LR particle.s; 
Col. 3: Mass of the DM component inside the virial radius; 
Col4:Numberofnearby clusters within 5 Hvir with M^M > 10^4^^-^ Mq. 



the adopted cosmology. Given the large volume this cosmological 
box contains a large sample of 64 clusters with Mfof > lO^'^/i"^ 
M0at 2 = 0. We selected the 24 most massive clusters for high 
resolution re-simulations. Figure lAll shows the projected density 
within 125 h^^ Mpc slices of the parent simulation at z=0. The po- 
sitions of the 24 most massive clusters used in this work are marked 
by diamonds. From the final output of the DM only run, all of the 
particles out to a distance of ~ 5 — 7i?vir around the cluster center 
were selected and then traced back to their initial positions. The 
corresponding Lagrangian region was enclosed in a box of side 
Lhr ~ 62.5 Mpc, the high resolution (HR) region. Since the vol- 
ume occupied by the HR particles, 1/hr, is usually only a fraction 
of the volume of the box (L^p), we sampled the box with 64"^ cells, 
and we marked cells which were actually occupied by the particles. 
In order to obtain a volume with a concave shape and no holes in 
it, some more cells were marked around/within Vhr- The particles 
that occupy the marked cells were then traced back to the initial 
redshift of the simulation. The right panel of Figure I A3 1 shows a 
cut through the Lhr volume. The blue cells trace the Vhk region, 
while the additional cells marked to obtain a concave volume are 
marked in red and green. This volume (defined as "occupied vol- 
ume") was re-sampled with a higher number of particles in order to 
obtain a higher mass resolution (in this case of 1 x lO^h"^ AIq{oi 
DM particles). Parti cles were d isplayed according to a glass-like 
particle distribution l lWhitell996) . The HR particles were perturbed 
using the same power spectrum of the parent simulation, keeping 
the same amplitudes and phases. New fluctuations at smaller spatial 
scales were added, since smaller frequencies are now sampled by 



the higher resolution particles. The amplitude of the fluctuations are 
given by the theoretical power spectrum P{k) of the parent simula- 
tion, but extended to higher k. To minimize any changes in the tidal 
forces acting onto the high resolution region, we created a buffer 
region around the HR region, and sampled it with the same mass 
resolution as the parent cosmological simulation. The remaining 
volume of the simulation was re-sampled at lower resolution ac- 
cording to the following procedure: the density and velocity fields 
of the LR particles were re-sampled onto a spherical grid having 
constant angular resolution d9. The size of each cell dr — rdO 
was chosen to obtain approximately cubic cells through the sphere. 
The interpolation onto a spherical grid reduces the number of LR 
particles to the minimum necessary to preserve the large-scale tidal 
field of the original simulation. We used d9 — 1.5°, resulting in 
~ 2 X 10*' low resolution par ticles, that guarante es an accurate 
sampling of the tidal field (see lTormen et al.ll 19970 . By construc- 
tion, as the distance from the HR region increases, dr increases 
too, and the mass of the LR particles increases accordingly. The 
overall volume simulated for each cluster is the same as the parent 
simulation, ensuring that the forming structures correspond to the 
same that formed within the original cosmological simulation. The 
new initial conditions were finally traced back to a higher redshift 
(e.g.z — 70) to ensure that the rms of the particle displacement 
in the HR region is still small enough to guarantee the validity of 
the Zeldovich approximation. After generating the new IC at higher 
resolution, we re-run further dark matter-only re-simulation to ver- 
ify that the volume of the HR region around each cluster was free 
from contamination of LR particles. Several iterations (typically 5- 
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7) of the whole procedure were reqmred for each cluster to obtain a 
clean, high resolution spherical volume with radius of 5 -Rvir , while 
keeping the total number of high resolution particles as low as pos- 
sible. In several cases, additional clusters close to our target had to 
be included in the high resolution region. Hence, all the initial con- 
ditions have at least a spherical volume of radius 5 — 6i?vir "clean" 
of low resolution particles, and centered on the target cluster (see 
Table lAlV The total number of high resolution particles needed is 
typically only 2 — 3 times larger than the number of high resolution 
particles within this regions of interest. Two of the selected clus- 
ters turned out to have a close-by companion with a mass larger 
than W^^h~^M(7). Whereas the 24 targeted clusters represent a 
fair volume-limited sample of galaxy clusters, the whole simula- 
tion sample encompasses in total 26 clusters with masses above 
lO^^h^^MQ. In addition, many other clusters with masses be- 
tween lO'^*/i~^A'/0and 10^^ h~^ Mq-wsk found close to our mas- 
sive targets. 50 of them are free from low-resolution particles up to 
at least 1 i?vir. We also extracted initial conditions of 5 more iso- 
lated cluster, having masses of « 5, 7, 4, 1, 1 x W^'^h~^ Mq. Such 
additional clusters are of interest when studying scaling relations 
( lFabianetani201lh . 



A2 Adding the baryonic component 

Once the IC for the DM particles have been obtained, the bary- 
onic component was added. The high resolution dark matter par- 
ticles are splitted into one gas and one DM particle. The mass of 
the initial DM particle is splitted according to the cosmic baryon 
fraction, conserving the center of mass and the momentum of the 
parent DM particle. We displaced them by half the mean interpar- 
ticle distance. Here we added a further optimization. Taking the 
"cleaned region" around all clusters of interest within the high reso- 
lution region, we traced back the corresponding Lagrangian region 
into the initial conditions. To associate concave volume to the se- 
lected particles, we measured their distance to the center of the high 
resolution region and calculated the maximum distance found in 
each direction by samp ling the sphere using a HealPIX discretiza- 
tion JGorski et alj2005n . Only those dark matter particles which are 
found within such a volume (including a very small safety buffer) 
were splitted into one gas and one DM particles. This typically 
saves « 20% of gas particles while still having splitted dark matter 
(and, accordingly, gas particles) within the full extent of the "clean 
region". 

In the left panel of figure IA3I in the central part, the spatial ex- 
tent of the whole high resolution region compared to the extent of 
baryon-filled region is visible. 
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Figure Al. 125 Mpc/i ^ thick slices through the parent DM simulation at 2: = showing the projected density. The Diamonds indicate the positions of the 
24 most massive clusters extracted for high resolution re-simulations. The bottom right panel shows the projected density through the whole box. 
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Figure A2. Ray-tracing images of a 15 Mpch ^ regions around the center of the individual clusters. Color coded is the temperature of the gas. 
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Table A2. Clusters with mass lQ^^h-'^MQ> M^^^ > lO^^h-^MQfree of low-resolution particles. 



Cluster name 


Mvi, 


Mgas 


i?vir 


^cleaned 




M0h-i 


h^'^MQ 


kpch-1 


Rvir 


dl.9 


1.436E-I-14 


2.067E-I-13 


1117.11 


>5 


d2.2 


1.708E-I-14 


2.506E-I-13 


1183.26 


>5 


d2.5 


1.646E-I-14 


2.331E-I-13 


1168.81 


>5 


d2.6 


1.070E-I-14 


1.540E-I-13 


1013.07 


5=5 


d3.4 


5.249E-I-14 


7.466E-I-13 


1723.11 


>5 


d3.23 


1.219E+14 


1.709E+13 


1057.81 


>5 


d5J 


7.7()7E-l-14 


l.l()4E-l-14 


1955.10 


4 


d5.6 


1.768E-I-14 


2.642E-I-13 


1197.06 


>5 


d5.ll 


2.266E-I-14 


3.327E+13 


1300.91 


4 


d5.25 


2.259E+14 


3.301E+13 


1299.61 


4 


d6.6 


1.434E-I-14 


2.117E+13 


1116.47 


>5 


d6.11 


2.286E+14 


3.191E+13 


1304.73 


>5 


d6.18 


1.310E+14 


1.971E+13 


1083.39 


1 


d6.26 


5.762E+14 


8.107E+13 


1777.52 


>5 


d7.11 


1.520E+14 


2.232E+13 


1138.18 


>5 


dS.l 


4.884E-I-14 


7.196E-I-13 


1682.05 


3 


d8.6 


4.993E-I-14 


6.883E-I-13 


1694.58 


>5 


d8.8 


2.112E-I-14 


2.993E+13 


1270.59 


2 


d8.29 


1.081E+14 


1.485E+13 


1016.29 


>5 


dlO.2 


1.838E+14 


2.532E+13 


1212.96 


>5 


dlO.3 


8.074E+14 


1.131E+14 


1984.74 


3 


dlO.4 


2.889E+14 


4.126E-I-13 


1411.13 


3 


dlO.6 


1.233E+14 


1.802E-I-13 


1061.86 


>5 


dl0.12 


1.017E-I-14 


1.475E-I-13 


996.10 


>5 


dlO.16 


1.057E+14 


1.346E-I-13 


1008.81 


>5 


dll.3 


1.119E+14 


1.562E+13 


1028.04 


>5 


dl2.1 


2.622E+14 


3.596E+13 


1366.08 


>5 


dl2.4 


3.815E+14 


5.462E+13 


1548.76 


>5 


dl3.1 


4.930E+14 


6.871E+13 


1687.93 


>5 


dl3.2 


3.808E+14 


5.518E+13 


1548.19 


>5 


dl3.3 


4.868E+14 


7.024E+13 


1680.78 


1 


dl3.7 


2.426E+14 


3.830E+13 


1331.19 


>5 


dl4.2 


1.344E+14 


2.043E+13 


1092.54 


>5 


dl4.3 


3.065E+14 


4.416E+13 


1439.40 


>5 


dl4.5 


1.754E+14 


2.519E+13 


1194.04 


^5 


dl5.7 


1.419E+14 


2.037E+13 


1112.67 


3:5 


dl7.4 


1.072E+14 


1.501E+13 


1013.47 


1 


dl8.1 


2.982E+14 


4.633E+13 


1426.32 


>5 


dl9.5 


2.573E+14 


3.267E+13 


1357.47 


>5 


d20.2 


6.489E+14 


8.920E+13 


1849.37 


>5 


d20.4 


1.135E+14 


1.615E+13 


1033.07 


>5 


d21.2 


2.083E+14 


2.978E+13 


1264.69 


4 


d21.3 


3.019E+14 


4.324E+13 


1432.09 


3 


d21.19 


1.031E-I-14 


1.687E-I-13 


1000.69 


>5 


d21.23 


1.526E-I-14 


2.141E-I-13 


1139.66 


>5 


d23.2 


1.134E+14 


1.619E+13 


1032.80 


>5 


d23.4 


2.970E+14 


4.217E+13 


1424.53 


>5 


d23.7 


1.172E+14 


1.681E+13 


1044.22 


5=5 


d24.1 


7.822E-I-14 


1.096E+14 


1964.50 


3 


d24.22 


1.149E-I-14 


1.485E-I-13 


1037.26 


>5 


d24.363 


3.341E+14 


4.041E+12 


1479.48 


>5 



Col. 1 : Cluster name; Col. 2: Total mass of the cluster inside the virial radius 
Col. 3: Mass of the gas component inside the virial radius; Col 4: Virial radius; 

Col 5: Number of virial radii cleaned by LR particles. 
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Figure A3. Initial condition region for the high resolution simulations. Left: Black: DM particles with degraded mass resolution outside the HR region, 
e.g.grained version of the original IC region used in the parent simulation, with increasing mass toward the outer regions. Green: DM particles outside the HR 
region with the same mass resolution than the parent simulation. This represents a "safety region" where a normal grid is used and particles have the same 
mass that the parent simulation. Red: HR region. Blue: region where high resolution DM particles have been splitted into gas and DM particles. Right: A slice 
through the HR initial condition region. Blue boxes refer to the position of the particles traced back, which where at 2 = falling within 5 i?vir of th^ target 
cluster. Red boxes are the cells that are included automatically to obtain a concave region. The green box refer to cell which was added by hand to avoid holes 
within the HR region. 



